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1. Introduction and General Status 

The main goals of the research under this grant consist of the development of 
mathematical tools and measurement of transport properties necessary for high fidelity modelling 
of crystal growth from the melt and solution, in particular for the Bridgman-Stockbarger growth 
of mercury cadmium telluride (MCT) and the solution growth of triglycine sulphate (TGS). Of 
the tasks described in detail in the original proposal, two remain to be worked on: 
development of a spectral code for moving boundary problems, 
diffusivity measurements on concentrated and supersaturated TGS solutions. 

During this eighth half-year period, good progress has been made on these tasks. 

2. MCT Code development 

During the last six monthly period we have completed a paper on a Chebyshev 
pseudospectral collocation method is adapted to the problem of directional solidification. 
Implementation of this method involves a solution algorithm that combines domain 
decomposition, a finite-difference preconditioned conjugate minimum residual (PCMR) method 
and a Picard type iterative scheme. The method solves equations which describe heat transfer in 
the ampoule, melt and crystal, and the convective flow problem in the melt. The crystal-melt 
interface shape is determined as part of the solution. A pre-print of the paper is attached in the 
appendix. 

In addition we have successfully completed an extension of this code to include species 
transport and the dependence of crystal melting temperature on composition. The method 
employs a conjugate-gradient-squared (CGS) technique for the species and heat transport 
equations, and a PCMR method for the momentum equations and involves finite-difference 
preconditioning. The code has been tested extensively against results of Kim and Brown [1] and 
Adomato and Brown [2] for the directional solidification of mercury cadmium telluride, gallium- 
doped germanium and silicon-germanium. Further work, beyond the tests, has involved the 
study of the interplay between convective flow, interface shape and compositional uniformity . 
These results will be reported in full in the next report. 

[1] D.H. Kim and R.A. Brown, “Models for convection and segregation in the growth of 
HgCdTe by the vertical Bridgman method”, J. Crystal Growth, 96, 609-627 (1989). 

[2] P. M. Adomato and R. A. Brown, “Convection and segregation in directional 
solidification of dilute and non-dilute binary alloys”, J. Crystal Growth, 80,155-190 
(1987). 
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3. Diffusivity Measurements 

Work on this topic has concentrated during the last six months on the accuracy of the 
novel diffusivity measurement technique developed under this grant. This was triggered by our 
discovery of poor reproducibility between runs. Two error sources were identified: 

The standard microscope slides used as windows in the diffusion cell, inspite of 
background interferogram subtraction, turned out to be optically inadequate to fully utilize the 
advantages of this technique. Hence, we have acquired optical windows flat to within 1/1 Oth of a 
wavelength of the He-Ne line used in the interferometry. This has led to a significant reduction 
of the experimental errors. 

In addition to the experimental error, we discovered that the mathematical approach taken 
in the evaluation of the interferometric data, can introduce larger errors than we expected earlier. 
Both the ZAPP-PC software used, as well as the specific function used in evaluating the integral 
equation used in our approach (see earlier reports) have been identified as significant error 
sources. Current work is concentrating on developing a more advantageous algorithm for data 
evaluation 

4. Presentations and Publications 

From the work carried out under this grant the following papers have been published, 
accepted for publication or are in preparation for submission for publication: 

1. A. Nadarajah, F. Rosenberger and J. I. D. Alexander, Modelling the Solution Growth of 
Triglycine Sulfate in Low Gravity, J. Crystal Growth 104 (1990) 218-232. 

2. F. Rosenberger, J. I. D. Alexander, A. Nadarajah and J. Ouazzani, Influence of Residual 
Gravity on Crystal Growth Processes, Microgravity Sci. Technol. 3 (1990) 162-164. 

3. J. P. Pulicani and J. Ouazzani, A Fourier -Chebyshev Pseudo-Spectral Method for Solving 
Steady 3-D Navier-Stokes and Heat Equations in Cylindrical Cavities, Computers and 
Fluids 20(1991) 93. 

4. J. P. Pulicani, S. Krukowski, J. I. D. Alexander, J. Ouazzani and F. Rosenberger, 
Convection in an Asymmetrically Heated Cylinder, Int. J. Heat Mass Transfer 35 (192) 
2119. 

5. F. Rosenberger, J. I. D. Alexander and W.-Q. Jin, Gravimetric Capillary' Method for 
Kinematic Viscosity’ Measurements, Rev. Sci. Instr. 63 (192) 269. 

6. A. Nadarajah, F. Rosenberger and T. Nyce, Interferometric Technique for Diffusivity 
Measurements in (Supersaturated) Solutions, J. Phys. Chem (submitted). 

7. F. Rosenberger, Boundary Layers in Crystal Growth, Facts and Fancy, in Lectures on 
Crystal Growth, ed. by H. Komatsu (in print). 
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8. F. Rosenberger, Short-duration Low-gravity Experiments - Time Scales, Challenges and 
Results, Microgravity Sci. Applic. (in print). 

9. Y. Zhang, J.I.D. Alexander and J. Ouazzani, A Chebishev Collocation Method For 
Moving Boundaries, Heat Transfer and Convection During Directional Solidification, 
Intemat. J. Numerical Methods Heat Fluid Flow (submitted). 

In addition to the above publications, the results of our work have been presented at the 

following conferences and institutions: 

1. J.I.D. Alexander, Modelling the Solution Growth of TGS Crystals in Low Gravity, 
Committee on Space Research (COSPAR) Plenary Meeting, The Hague, Netherlands, 
June 26- July 6, 1990. 

2. A. Nadarajah, Modelling the Solution Growth of TGS Crystals in Low Gravity, Eighth 
American Conference on Crystal Growth, Vail, Colorado, July 15-21, 1990. 

3. J.I.D. Alexander, Commercial Numerical Codes: To Use or Not to Use, Is This The 
Question? , Microgravity Fluids Workshop, Westlake Holiday Inn, Cleveland Ohio, 
August 7-9, 1990. 

4. F. Rosenberger, Fluid Transport in Materials Processin, Microgravity Fluids Workshop, 
Westlake Holiday Inn, Cleveland Ohio, August 7-9, 1990. 

5. F. Rosenberger, Influence of Residual Gravity on Crystal Growth Processes, First 
International Microgravity Congress, Bremen, September 1990 (invited). 

6. J.I.D. Alexander, Residual Acceleration Effects on Low Gravity Experiments, Institute de 
Mechaniques des Fluides de Marseilles, Universite de Aix-Marseille III, Marseille, 
France, January 1991, (3-lecture series, invited). 

7. J.I.D. Alexander, An Analysis of the Low Gravity Sensitivity of the Bridgman-Stockbarger 
Technique, Department of Mechanical Engineering at Clarkson University, April 1991 
(invited). 

8. A. Nadarajah, Measuring Diffusion Coefficients of Concentrated Solutions, Fifth Annual 
Alabama Materials Research Conference, Birmingham September 1991. 

9. A. Nadarajah, Modelling Crystal Growth Under Low Gravity, Annual Technical Meeting 
of the Society of Engineering Science, Gainesville, November 1991. 

10. J.I.D. Alexander, Vibrational Convection and Transport Under Low Gravity Conditions, 
Society of Engineering Science 28th Annual Technical Meeting, Gainesville, Florida, 
November 6-7, 1991. 
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11. F. Rosenberger, Theoretical Review of Crystal Growth in Space - Motivation and Results, 
International Symposium on High Tech Materials, Nagoya, Japan, November 6-9, 1991 
(plenary lecture, invited). 

12. F. Rosenberger, Computer Simulation in Materials Science, Mitsubishi Frontiers 
Research Institute, Tokyo, Japan, November 8, 1991 (invited). 

13. F. Rosenberger, Importance of Materials Research in Space Laboratories for Industrial 
Development, International Symposium for Promoting Applications and Capabilities of 
the Space Environment, Tokyo, Japan, November 14-15, 1991 (plenary lecture, invited). 

14. F. Rosenberger, What Can One Learn from 10 Second Low-Gravity Experiments?, In 
Space 1991, Tokyo, Japan, November 14-15, 1991 (plenary lecture invited). 

15. P. Larroude, J. Ouazzani and J.I.D. Alexander, Flow Transitions in a 2D Directional 
Solidification Model, 6th Materials Science Symposium, European Space Agency, 
Brussels, Belgium, 1992 (poster). 

16. F. Rosenberger, Microgravity Materials Processing and Fluid Transport, AIAA Course 
on Low-Gravity Fluid Dynamics, AIAA Meeting, Reno, NV, January 10-12, 1992 (3- 
lecture series, invited). 

17. J.I.D Alexander, Numerical Simulation of Low-g Fluid Transport, AIAA Course on Low- 
Gravity Fluid Mechanics, Reno, NV, January 10-12, 1992 (invited). 

18. F. Rosenberger, Time Scales in Transport Processes and Challenges for Short-Duration 
Low-Gravity Experiments, Falltower Days Bremen, Bremen, Germany, June 1-3, 1992 
(invited). 

19. J.I.D. Alexander, Modelling or Muddling? Analysis of Buoyancy Effects on Transport 
under Low Gravity Conditions, World Space Congress, Washington, DC, August 28 - 
September 5, 1992 (invited lecture). 


A CHEBYSHEV COLLOCATION METHOD FOR MOVING 
BOUNDARIES, HEAT TRANSFER, AND CONVECTION 
DURING DIRECTIONAL SOLIDIFICATION 


Yiqiang Zhang, J. Iwan D. Alexander and Jalil Ouazzani* 

Center for Micogravity and Materials Research, University of Alabama in Huntsville 

Abstract 

Free and Moving Boundary problems require the simultaneous solution of unknown field variables 
and the boundaries of the domains on which these variables are defined. There are many 
technologically important processes that lead to moving boundary problems associated with fluid 
surfaces and solid-fluid boundaries. These include crystal growth, metal alloy and glass 
solidification, melting and flame propagation. The directional solidification of semi-conductor 
crystals by the Bridgman-Stockbarger method 1 - 2 is a typical example of a such a complex process. 
A numerical model of this growth method must solve the appropriate heat mass and momentum 
transfer equations and detemine the location of the melt-solid interface. In this work, a Chebyshev 
pseudospectral collocation method is adapted to the problem of directional solidification. 
Implementation of the method involves a solution algorithm that combines domain decomposition, 
a finite-difference preconditioned conjugate minimum residual method and a Picard type iterative 
scheme. 


* Presently at the Institute de Mceanique des Fluides de Marseille, 1 rue Honnorat, Marseille, France. 
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1. INTRODUCTION 

Moving and free boundary problems are problems that require as part of the solution the 
determination of some or all the boundaries of the domain under consideration. Included in this 
class of problems are situations that involve fluid surfaces, or solid-fluid interfaces. Freezing and 
melting, crystal growth, flame propagation, liquid surface configurations, are examples of such 
processes that are important in a variety of areas with technological applications. Such problems 
generally pose a challenging problem to the numerical modeller. The Bridgman-Stockbarger 
directional solidification crystal growth technique is a typical example of such a complex problem. 
To adequately represent the physics of the problem, the solution method must be able to cope with 
the following: The unknown location of the crystal-melt interface, high Rayleigh number 
buoyancy-driven flows, heat transfer by conduction (along ampoule walls and in the crystal), 
convective-diffusive heat transfer in the melt and radiative and convective heat transfer between the 
furnace and the ampoule. Even for pure melts, due to differences in thermal conductivities between 
melt , crystal and ampoule, and the differences in thermal and momentum diffusivities in the melt, 
the problem has a variety of disparate length scales over which characteristic features must be 
accurately represented. 

In past work 3 ' 11 , the Finite Element Method (FEM) has been successfully applied to the 
problem of computing melt and crystal temperature and concentration distributions, melt 
convection and the location of the crystal-melt interface. As an alternative to FEM we present a 
Chebyshev collocation (pseudospectral) method suitable for the solution of this class of problem. 
Spectral and pseudospectral methods *2- 13 involve the representation of the solution as a truncated 
series of smooth functions of the independent variables. In contrast to FEM, for which the 
solution is approximated locally with expansions of local basis functions, spectral methods 
represent the solution as an expansion in global functions. In this sense they may be viewed as an 
extension of the separation of variables technique applied to complicated problems 14 . 

For problems that are characterized either by irregularly shaped domains, or even domains 
of unknown shape, it is, in general, neither efficient nor advantageous to try to find special sets of 
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spectral functions that are tuned to the particular geometry in consideration (especially in the case of 
solidification, where the melt-crystal geometry is not known a priori ). Two alternative methods 
are mapping and patching 14 . Mapping allows an irregular region to be mapped into a regular one 
(which facilitates the use of known spectral functions, such as Chebyshev polynomials). For 
directional solidification systems (see Fig. 1) the melt-crystal boundary and, thus, the melt and 
crystal geometries, are unknown. Nevertheless, by specifying the melt-crystal boundary as some 
unknown single-valued function, the melt and crystal geometries can be mapped into simple ones 
by a smooth transformation. This mapping facilitates the use of Chebyshev polynomials to 
approximate the dependent variables in these new domains. 

As can be seen from Figs. 1 and 2, heat transfer to and in the ampoule wall must also be 
considered. To do this we employ patching by subdividing the system into four domains (crystal, 
melt and two ampoule domains), and transform these domains to domains with simple shapes. We 
then solve the resulting problems in each domain and solve the full problem in the complicated 
domain by applying suitable continuity conditions across any boundaries (real or artificial) between 
the domains. 

The formulation of the problem is outlined in section 2. The solution method is described in 
section 3. Our results are presented in section 4 and discussed in section 5. 

2. FORMULATION 

The vertical Bridgman-Stockbarger system is depicted in Fig. 1. A cylindrical ampoule 
with inner and outer diameters of 2Ro and 2(Ro+R w ) contains melt and crystal. To grow the crystal 
the ampoule must be translated relative to a prescribed external temperature gradient. The objective 
of this model is to describe a steady growth process that, in reality, can be achieved between initial 
and terminal transients in sufficiently long ampoules. Toward this end a pseudo-steady state 
model 2 is employed that neglects the ends of the ampoule. The remainder of the ampoule is 
assumed to occupy a cylindrical computational region of length L. Ampoule translation is then 
accounted for by supplying a melt to the top of the computational space at a uniform velocity, and 
withdrawing crystal from the bottom at the same velocity. It is thus assumed that there is no lag 
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between the translation rate and the crystal’s growth velocity. Transport of heat from the furnace to 
the ampoule is modelled using a prescribed furnace temperature profile. The heat transfer from the 
furnace to the outer ampoule wall is governed by a heat transfer coefficient Bi(z). This is discussed 
in more detail later. The top and bottom of the ampoule are respectively assigned temperatures of 
Tr and Tc (Th > Tc). 

The variables are cast in dimensionless form by using Ro, cci/Ro ctTRo 2 anc * 

Th-Tc, where ocl is the melt’s thermal diffusivity, to scale length, velocity, stream function, 

vorticity and temperature, respectively. That is, 

? T T 

x = (r,z)=( r,z)/Ro u = u R</a L , \\f =V/ Root , © = ®Ro/a L , T = - ' S - ■ (1) 

Th-Tc 

Here r and z represent the radial and axial coordinates, \\r is the stream function, © is the vorticity 
and u = (u,w) represents velocity with radial and axial components u and w, respectively. A tilde 
denotes a dimensional quantity. Melt, crystal and ampoule temperatures will be distinguished by 
the suffixes L (melt), S (crystal)and W (ampoule) when necessary. The location of the crystal melt 
boundary is given by z=h(r,t) and must be determined. The melt is assumed to be incompressible 
and the stream function and vorticity are defined by the velocity components (u,w) as 


_ 1 

U r dz 


w 


__1 d\\f _9u _ dw 
"Tl)? 03- ^ dr ’ 


The governing equations then take the following form 
In the melt, 0 <r < 1 , 0< z < h(r,z) 


( 2 ) 


u ^ +w ^-_miI=Pr 

dr dz r 


r 8 2 co ] 3© 3 2 © 

^ + dz 2 


- Pr — + Pr 2 Gr^ 
r 2 3r 


rw = 


d\r i dy_ 3V 
~d?~T dr + 3z 2 ’ 


( 3 ) 

( 4 ) 


and 
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ai ii.ii 

dr dz ’ dr 2 r dr dz 2 


dT dT 
u — + w— = AT 


( 5 ) 


where Pr = v/ocl is the Prandtl number, Gr = (3(T H - Tc)gRo 3 /v 2 is the Grashof number, v is the 
melt’s viscosity and (3 is the melt’s thermal expansion coefficient. 


In the crystal, 0 < r < 1 , h(r) « z < A , 

D I FT. 

a Pe^r — = AT, 
dz 


( 6 ) 


and in the ampoule wall, 1 < r < r w , 0< z < A, 


a'Pe^— = AT, (7) 

dz 

where a', and a" are, respectively, the ratios of the melt’s thermal diffusivity with the crystal and 
ampoule thermal diffusivities, and Pe = VqRo/06l is the Peclet number and Vo is the ampoule 
translation rate. 

For the temperature the boundary conditions are: 

At the melt-crystal interface z= h(r,t) 

T l = Ts = T m , (8) 

k'VT L n- V T s n=StPe a'n-e z , (9) 

where Tm represents the dimensionless melting temperature, k' is the ratio of melt and crystal 
conductivities, St = AH/(C ps AT) is the Stefan number. The vector n is the unit normal to the 
crystal-melt surface and points into the melt. At the outer ampoule wall, r= r w 

-^ = Bi(z)(T-T F (z)). (10) 

The temperatures at z=0 and z=A are constant, i.e. 


T(r,0)= 1, T(r,A) = 0, 
and the heat flux is continuous across the inner ampoule wall 


( 11 ) 
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/dT(l,z)\ 

= k JST(U)\ 

l dr ) L 

l dr I- 

[dT(l,z)\ 

.»,OT(uy 

1 dr )s 

l dr i 


( 12 ) 


w 


In (10) Bi(z) is a heat transfer coeffcient and Tp(z) is the furnace temperature profile. The 
coefficients k* and k** represent the ratio of the wall conductivity with that of the melt and 
ampoule, respectively. 

For the stream function the boundary conditions are 


\|/ (0,z) = 0, \|/(l,z) = - \ Pe, y(0,z) = - 1 r 2 Pe, \|/(h(r),z) = - ^ r 2 Pe, 


(13) 


and the vorticity is zero at r=0. At the other melt boundaries the boundary conditions for the 
vorticity are enforced (iteratively) using previously computed values of the velocity field (the 
scheme is explained in section 3.3. The velocity boundary conditions are 

u(0, z) = u( 1 ,z)=u(r,0) = u(r,h(r))=0 = ^ W ^’ Z) ; w(r,0) = w(l,z)= w(r,h(r))=Pe. (14) 

Note that, at the melt-crystal boundary there are two boundary conditions for the 
temperature. In the following section we describe an iterative scheme which distinguishes one of 
the temperature boundary conditions and uses it to compute the interface shape iteratively. 


3. SOLUTION METHOD 

The solution method is based on a Picard 15 type iteration which consists essentially of four 

steps: 

1. The initial shape of the crystal-melt interface is specified and an independent variable 
transformation is applied to the governing equations and boundary conditions in the melt, crystal 
an ampoule regions. This specifies the computational domains. 

2. The coupled momentum, heat, mass and species equations are then solved using three of the 
four boundary conditions on the moving boundary. 
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3. The remaining boundary condition (or distinguished condition 2 ), in this case equation (8), is 
used to compute corrected boundary locations. 

4. Steps 2 and 3 are repeated until the distinguished boundary condition is satisfied. 

The solution method is implemented using domain decomposition and a preconditioned 
generalized conjugate residual method 13 * 1 ** 

3.1 Domain decomposition 

The physical region is split into four computational domains, Qj, i=l,...,4. The domains 
correspond to the melt (Hi), the crystal (E3), and the portions of the ampoule wall adjacent to the 
melt (H2) and the crystal (E4), as shown in Fig. 2. The irregularly shaped domains are mapped 
onto rectangular regions by 


4=r ' ,1= h(V i '^ 

(15) 

s =r -1=h(l)' a ^ 

(16) 

^= r .11- 2 -h(r)-A.-3^ Q 3- 

(17) 

^ =r * T ' — 2 — h(l)— A’ 111, 4_> ' 

(18) 


3.2 Spatial discretization 

The dependent variables, <E> are approximated by Chebyshev polynomials 12 * 13 , i.e 

N,M 

<t>=<t NM (X i ,Y J )=X a ijTijlXj.Vj), (19) 

i = 0 
j = 0 

where Tjj = TjTj , and the Tk are Chebyshev polynomials of order k. The points (Xj , Yj ) are 

related to the coordinates c, and r) by 

^ = aX + b, Tj = cY + d, (20) 

where a and b are determined by the transformation of each domain, £2j, to [-1,1 ]x[— 1 , 1 ]. The 
discrete points (Xi,Yj), i=0, N, j=0,M,are the Gauss-Lobatto collocation points 13 . That is, 
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X : = COS K 


X; = COS 71 


1 

N 

i 

N 


,i = 0. 


,i = 0, 


The spatial derivatives are given by 

ao i ao ao i ao a 2 o _ j a 2 o a 2 <J>_ j a 2 o> a 2 o_ i a 2 o 
a^ " a ax’ a^ " c aY ’ 3PT “ ac axay’ a^ 2 a 2 ax 2 ’ a^i 2 c 2 av 2 

where the derivatives with respect to X and Y have the forms 

£ Dj^Yj). £ DiPOpj, 

dX 3X p = o p = o 

ao(x, v ) ,9Q(x „ v j j L= £ 0>(Xi , Yq|= £ # iq , 

aY aY q = 0 q = 0 


ax" 


ax p = o 


p = 0 


a 2 o(x,Y)_a 2 a)NM(x i ,Y j ),_ 


dY Z 


M M 

£ D*4>(Xi,Y,)=£ D j? y <K iq. 

aY q = 0 


q = 0 


N, M 


N, M 


3<1> < Xi ’ ^ = £ D i M‘l<I>(X p ,Y q ) = £ 
aXY p = o P = o 

q = 0 q = 0 


( 21 ) 

( 22 ) 

(23) 


(24) 


(25) 


fr(X,Y) .i»<*Y,X. £ D iP 0(Xp , Yj)= £ D .PO pj, (26) 


(27) 


(28) 


where the expressions for D x , D y , D xx , D yy and D xy are given explicitly by Ouazzani 17 . 

3.3 Pseudo-unsteady discretization 
The governing equations now have the form 

(A«-A (i) >l) (i) = S« i= 1,6 (29) 

where the A(»), A( i ), S(') and <)>(•) are given in the Appendix B. To solve these equations using a 
pseoudo- unsteady iterative scheme, equation (29) is rewritten as 
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(^ + A«}(t) (i) ==S'' i ) + A (i) (t) (i) , i= 1,6 

ol 


( 30 ) 


The left-hand side of (30) is written in discrete form as 

(^+ A )<[> =(A + al)<|) n+1 — cr<t> n , 
or 


(31) 


where the index in parentheses has been omitted for clarity and the superscript denotes the pseudo- 
time or iterative step number. Note that, cKO = 1/Ax(‘) for i=l-6 and that the time step size, At( 0, is 
generally different for each of the equations. 

The problem can now be expressed as 

H sp <!> n+1 + al (|) n+1 = F(<t>,h) n , (32) 

where 

H sp = A -A, and F= S + a<|> n , (33) 

and H sp j s obtained from the expressions in Appendix B using the Chebyshev derivatives (24)- 
(28) and equation (23). A superscript n denotes a quantity evaluated at the nth iterative step (note 
that the indices in parethenses have been omitted for clarity) . 

3.4 Vorticity boundary condition 

To solve the vorticity-stream function equations we adopted the following procedure which 
is simply an extension, for Chebyshev approximations, of an approach described by Peyret 1 -. The 
velocity field is calculated from the stream function obtained from the previous iteration. The 
vorticity at the boundary which corresponds to this velocity field is then found from 


-n+l_pu _8w\ n 
\dz dr r 


(34) 


and the value of the vorticity to be applied at the boundary, co n+l , is given by 

to n +] = Y» n+1 +(l-y)0) n . 


(35) 


Here y(0 < y< 1) is a relaxation parameter. 
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3.5 Preconditioned generalized conjugate residual method 

The operator H sp is represented by a full matrix of order (N+l) 2 x (M+l) 2 and is not 
symmetric. In order to solve the system of equations and boundary conditions represented by (32) 
and (A.l 1 )-(A. 19), each of the spectral operators H sp for each of the domains Hj , i= 1 ,4 and the 
conditions on shared domain boundaries Qj n Qj b*j, i,j=l,4 are combined and approximated by a 
single finite difference operator Hfy. This is defined over the entire domain The 

following iterative procedure which consisits of inner and outer loops is then adopted: 

Outer loop: First an initial interface shape h° is assumed 
Inner loop: The residual R is then initialized by 



where O represents the <f>W. Then we solve 

H; d O° = R°, 


where H = H + ol. Then we set 


P° = 0°, 


and calculate 


^m+l 


(R m ,H sp P m ) 

*rn it nrrK 


(H sp P >H sp P ) 

The variables O are then updated from 

<D m+1 =«D m + a m+1 P m , 


and the problem 


H; d © m+, =R m+i , 


(36) 


(37) 


(38) 


(39) 


(40) 


(41) 


is solved for 0. P is then updated using 
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p mf1 =0m+ £ pjwlpl, 

J = 0 J 


( 42 ) 


where 

PfH^— sp — : 2_. (43) 

<H*Ptty“) 

The procedure is continued until | R | < E. 

The preconditioned problem is given by equations (37) and (41). The finite difference 
operator H*fy is approximated by incomplete LU decomposition. The solution for © is obtained 
by forward and backward substitution. The subsequent approximations to d> = (Ts,Tl,cd,\|/ ) are 
then obtained from (40). At this point we note that while we used a nine-diagonal matrix for the 
second-order central finite difference operator for the solution of the temperature field, a seven 
diagonal operator was used for the solution of the stream-function and vorticity as it appeared to 
lead to more rapid convergence. This means that the cross-derivative terms were evaluated at the 
previous time step and were included in F on the right-hand side of (32). 

3.6 Interface shape update 

This iterative procedure is repeated until the convergence criterion is satsified. The first of 
equations (A. 17) is used as a distinguished boundary condition. If it is not satisfied, another outer 
loop iteration is performed and the interface shape is relocated using either Newton’s method 

h " +1 =h' 1 + ||^j" 1 e i (44) 

where 6j is the difference between the temperature at the ith interfacial site and the melting 
temperature T m ; or from a searching method 


hf +l =h" + a(Tj-T M ). 


(45) 
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Here a is found by numerical experiment. We found that by using the Newton method for the first 
few iterations and then the searching method for subsequent iterations, we achieved better success 
than with the Newton method alone. 

4. RESULTS 

We carried out several tests of the method. The results are shown in Tables 1-3 and in 
Figs. 3 and 4. The parameters used are given in Appendix C and correspond to the thermophysical 
properties of Gallium-doped Germanium. For the cases examined our results are in good 
agreement with the FEM calculations of Adomato and Brown.*' (see table 2). 

Figure 2 shows results for a furnace with a constant temperature gradient and Bi=7.143. 

That is, 

Tf(z) = 1 -zA 1 . (46) 


The isotherms are practically flat except at the crystal-melt boundary where the mismatch in thermal 
conductivity results in a convex interface. The flow depicted by the streamlines in Fig. 3b-d is a 
combination of the ampoule translation (which, if buoyant convection were absent, would appear 
as a set of vertical streamlines parallel to the ampoule wall) and buoyant flow caused by radial 
gradients in temperature. This results in an downward flow of hot melt near the ampoule wall and a 
upflow near the ampoule centerline. Note the increase in flow intensity as the Grashof number is 
increased. 

Figure 4 shows results for different Grashof numbers for a non-uniform furnace 
temperature profile given by 

T((z) =0.5[1+ tanh (6-12zA~’)] (47) 

together with a position dependent heat transfer function given by 

Bi(z) =0.2{2[1+ tanh (5-2z)] + 1+ tanh (2z-15)}. (48) 

Radial temperature gradients arise for two reasons in this problem: The mismatch in thermal 
conductivities at the ampoule-melt-crystal junction and the change in heat transfer at the quasi- 
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adiabatic zones. These zones are created by the furnace temperature profile and conditions (47) and 
(48). This heating configuration produces two counter rotating cells. The upper cell increases in 
spatial extent as the Grashof number is increased. 

Table 1 shows the CPU times, number of iterations taken to converge and compiler options 
for the case shown in Fig. 4b for a CRAY/XMP, an iPSC parallel processor and an Ardent Titan 
computer. 

5. DISCUSSION 

Chebyshev spectral methods that have been shown to achieve superior accuracy for a wide 
range of fluid flow problems defined in regular geometries can be applied to problems involving 
unknown free and moving irregular boundaries through a combination of mapping and domain 
decomposition. For the directional solidification described here, this was achieved without 
incurring excessive CPU times and has been implemented on several different machines to 
illustrate the magnitude of the CPU times involved for a typical calculation. Whether there is 
ultimately any advantage in using such spectral methods over finite elements will depend on the 
specific application. It will most likely depend on the accuracy required and on whether the ability 
of the Chebyshev collocation method to achieve better accuracy for a given number of collocation 
points (which is recognized for a variety of flows in regular geometries) is retained or degraded 
when using domain decomposition. 
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Appendix A 
Transformed Equations 

After the equations and boundary conditions (2) - (14) have been transformed according to 
(15)-(18) we obtain the following equations. 

For 0 < r| < 1 
0 < E, < 1 


where 


and 


For 1 < T| < 2 
0 < £, < 1 


1 5T 5v|/ 5T 5vj/ . _ A *y 
& 5n 5n ar 

(A.l) 

j_ 3\|/3co _5\[/5co _co5t|/ p rA * r ^ _ p^ + PrGij^- 

%5n ^dn % w h drdr ^ 

), (A.2) 

3V +a 3V + b |X-I^ + c|M<o. 

3 ^ 2 5r| 2 5<;5ri \ 5s 5r| 

(A. 3) 

A*= 32 + A 32 + B 32 tflr+cA, 
9 ^ 2 5r| 2 ^ 5q 5r| 

(A.4) 

A 1 +H dhf B " 2T ldh C -n[2( ldh ) 2 l d2 h--Ldh.j 
A 'S lh dFr h dr ^HhdrJ h ifi dr j 

(A. 5) 


A ** T _ p e a ' 1±L_ = o , 

hdr\ 


(A. 6) 


where 


and 


A ** = iL +A *— + f Jr + c*J-, 

3^ 2 5q 2 5^5q £5/; 5r| 


A* = • 


(h-A) 


+ , JlzL dh r B * = ^2_ _2) dh , 

7 'h-A dr/ h-A dr 


C* =01 -2) 


4 1 dh\ 2 L 

Ih-A dr) h- 


d^i 


dh 


A dr 2 ^(h-A)dr 


(A. 7) 


(A. 8) 
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In the ampoule wall, 1 < £ < r w , 0 < p < 2, where h(r) is taken to be a constant at each inner 
iteration, we have 


and 


^L + i|L + l|-Pea|| = 0,«<n<l, 
a^ 2 h 2 ^n 2 £ ^ h 

— + — *— tl? ' +Pea 'iT L T^ : = 0 - 1<T1<2 

dt? (h-A^dp 2 h-Adn 


The boundary conditions become 


aT 


= 0, \|/ = 0, co =0 at ^ = 0, 


T = 1, \|/ = - atr) = 0, 


\|/ = -^-Pe, 


/dT L 

T] dh dT L \ 

U 

h dr dp | 


* 


/dT s 

P dh 3T s| 

\ d^ 

h - A dr dp ) 




^ = Bi(p)(T-T F (p)), ^ = r w ,0<p<2 , 

d^ 

T= 0 at t| = 2, 0 < \ < r w . 

Finally, at the crystal melt interface the boundary conditions are 


t-t 1 k^Tl . 1 3T S _ StPecx' 

M ’ l/dp h - A dp 1+ ^dh) 2 ’ 


(A. 9) 
(A. 10) 
(A. 11) 
(A. 12) 

(A. 13) 

(A. 14) 

(A. 15) 
(A. 16) 

(A. 17) 


and 
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\|/ = -^; 2 Pe, atT| = l,0<q< 1 


(A. 18) 


In (A. 17) we have used the fact that the melting temperature Tm is assumed to be constant along 
the crystal melt interface (i.e. dT/3^ = 0). The vorticity boundary condition is given by equation 

(35) with 

. 3 .,, 

(A. 19) 


_ 1 ■ *1 dh d w 

h dp h dr dr) d \ 


Appendix B 

The A(b, AO) and F0) referred to in section 3.3 are expressed in terms of the equations 

given in Appendix I as follows: 

<J)(I) = T n+1 


<J>(2) = to n+1 


(|)(3) = y 


rn+1 


Ad) _ 1 pr n+l dy n 

9i an ^ 

A (,) = A* 

Fd) = tfdjT n 


r2 i _ l 3co n+I avj/ n 3co n+1 (0 n+idi|/ n . 

an ' l 3n ’ 

a < 2 > = A* - Et 


F® = PiW^- - 5- f- 5^1 + ^2)co" 

\ h dr an I 


a (3) =0 

A (3) =A *_2i 
F (3 ) = ^(0 n+1 + a (3) \|/ n , 

((>(4) = T w n+1 (0 < p < 1) , 4> f5 ) = T w n+1 (1 < p < 2), = T s n+1 , 


(B.l) 

(B.2) 

(B.3) 

(B.4) 

(B.5) 

(B.6; 

(B.7) 

(B.8) 

(B.9) 
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A< 4 > =A® = A“> =0 

(B.10) 

A«> + 1 f + 1 — -Pea'tJ- 

£2 h 2a n 2 !;3£, h3ri 

(B.11; 

A< 5 > _^-+ 1 32 + 1 3 +Pea” 1 3 

3^ 2 (h-A) 2 ^ 2 h - A 

(B.12; 

A (6) = A**- Pe a' 7-=^- , 
hdri 

(B.13) 

F« = c(»<l> (i) , i = 4,5,6 

(B.14) 
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Appendix C 

Physical constants, system dimensions and thermophysical properties of Gallium doped 
Germanium used in the calculations 


Property 

dimension 

Ge:Ga 

Growth velocity 

[cms' l i 

4X10- 4 

Ampoule length (L) 

[cm] 


Constant gradient furnace (Fig. 2) 


7.0 

Heat pipe furnace (Fig. 3) 


7.62 

Outer ampoule radius (R w ) 

[cm] 


Constant gradient furnace 


0.7 

Heat pipe furnace 


0.952 

Inner ampoule radius (Ro) 



Constant gradient furnace 


0.5 

Heat pipe furnace 


0.762 

Kinematic viscosity (v) 

[cm 2 s'* ] 

1.3(10)‘ 3 

Thermal conductivity (ampoule) 

[WK-'cnr 1 ] 


Constant gradient furnace 


3.27 

Heat pipe furnace 


0.26 

Thermal conductivity (crystal) 

[WK-'cm- 1 ] 

0.17 

Thermal conductivity (melt) 

[WK^cnr 1 ] 

0.39 

Density (crystal) 

[g cm- 3 ] 

5.5 

Density (melt) 

[g cm- 3 ] 

5.5 

Heat of solidification (AH) 

[J S' 1 ] 

460 

Specific heat (melt) 

JK-lg- 1 

0.39 

Specific heat (crystal) 

JK-lg-l 

0.39 

Thermal expansion coefficient 

[K-'] 

5 (10)- 4 
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Figure Captions 

Fig. 1 Typical Bridgman-Stockbarger set-up 

Fig. 2 a) The model Bridgman-Stockbarger system and b) the computational domains 
Fig.3 Results for results for a furnace with a constant temperature gradient, Bi=7. 143 and a Pr = 
0.07 melt, a) Gr = 5206, b) Gr = 52,060 c) Gr = 520,600 
Fig. 4 Results for a non-uniform furnace temperature profile (47) and position dependent heat 
transfer coefficient (48) for Pr =0.007 and a) Gr = 7,140, b) Gr = 14280 c) Gr=7 1,400 
d) Gr= 142,800 . 












